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The ionization of atoms by strong, low-frequency fields can generally be described well by as¬ 
suming that the photoelectron is, after the ionization step, completely at the mercy of the laser 
field. However, certain phenomena, like the recent discovery of low-energy structures in the long- 
wavelength regime, require the inclusion of the Coulomb interaction with the ion once the electron 
is in the continuum. We explore the first-principles inclusion of this interaction, known as analytical 
i?-matrix theory, and its consequences on the corresponding quantum orbits. We show that the 
trajectory must have an imaginary component, and that this causes branch cuts in the complex 
time plane when the real trajectory revisits the neighbourhood of the ionic core. We provide a 
framework for consistently navigating these branch cuts based on closest-approach times, which 
satisfy the equation r(t) • v(t) = 0 in the complex plane. We explore the geometry of these roots 
and describe the geometrical structures underlying the emergence of LES in both the classical and 
quantum domains. 


The interaction of atoms and molecules with intense 
lasers is a rich field, and provides challenges for theory to 
describe non-perturbative phenomena which often bridge a 
large range of energy scales. A recent example of this is the 
discovery of low-energy structures (LES) [1, 2] in photoion¬ 
ization by strong, long-wavelength fields where, in addition 
to above-threshold electrons that absorb many more pho¬ 
tons than required to reach the continuum, a strong peak is 
observed at energies far below the mean oscillation energy 
of electrons in the field. 

This was unexpected [3], as the strong-field approxima¬ 
tion (SFA) [4] typically predicts a smooth, featureless spec¬ 
trum at low energies, and SFA is generally expected to im¬ 
prove in accuracy as the frequency decreases and the sys¬ 
tem goes deeper into the optical tunnelling regime. This 
triggered active interest in developing theoretical methods 
to describe the LES, and in identifying the physical mech¬ 
anisms that create them. 

Numerical simulations of the time-dependent Schro- 
dinger equation (TDSE) do reproduce the structure [1], 
though they are particularly demanding at long wave¬ 
lengths. Simulations done with and without the ion’s long- 
range Coulomb potential [5, 6] indicate that its role is es¬ 
sential in producing the LES. Similarly, Monte Carlo simu¬ 
lations [6-9] using classical trajectories involving both the 
laser and the Coulomb field support this conclusion. 

From a semiclassical perspective, it is indeed possible to 
include the effect of the ion’s potential when the electron 
is already in the continuum. This can be done via a per¬ 
turbative Born series, like those used to explain high-order 
above-threshold ionization [10]; here the LES emerges as 
electrons that forward-scatter once at the ion [11-14]. Al¬ 
ternatively, the SFA can be reformulated to include the ef¬ 
fect of the Coulomb field on the underlying trajectories [15- 
17]; the resulting Coulomb-corrected SFA (CCSFA) also 
points to forward-scattered electrons [18]. First-principles 
analytical methods to include the Coulomb field’s effect on 
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the wavefunction’s phase, known as analytical i?-matrix 
theory (ARM) [19-24], are as yet untested in this regime. 

Analysis of the classical trajectories involved in the scat¬ 
tering points more specifically at soft recollisions, where the 
electron does not hit the core head-on, but is instead softly 
deflected as it approaches the ion with nearly zero velocity, 
near a turning point of the laser-driven trajectory [25, 26]. 
Although these trajectories spend enough time near the 
ion that the effect of its potential is no longer perturbative 
(and can in fact cause chaotic dynamics [27]), a comparison 
of the different models suggests that the LES are rooted in 
the pure laser-driven orbits of the simple-man’s model, and 
that the role of the Coulomb potential is to enhance their 
contribution [28]. 


ZAqijjv 



FIG. 1. Trajectories with soft recollisions after tunnel ioniza¬ 
tion, for a Keldysh parameter of 7 = 0.75. 

This raises an intriguing possibility, because there is in 
general a discrete sequence of soft-recolliding trajectories 
[26, 28], shown in Fig. 1, which spend increasing amounts of 
time in the continuum before recolliding, and which should 
appear as distinct peaks in photoelectron spectra. In this 
work we extend this sequence and show that there are, 
in fact, two families of soft-recolliding trajectories, which 
revisit the ion at even or odd numbers of half-periods after 
ionization. 
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The previously-described trajectories [26,28-30] recol¬ 
lide after an odd number of half-periods, and their drift 
energy is a constant multiple of the ponderomotive poten¬ 
tial Up. In contrast, trajectories with soft recollisions at 
integral multiples of the laser period have lower energies, 
of around 1 meV for common experimental parameters, 
which scale rather unfavourably as 1 /U p . This series seems 
to have been overlooked, but the low energies mean that the 
peaks can (and, in fact, will) contribute to the ‘zero-energy 
structure’ (ZES) found recently [31-33]. 

Moreover, we show that both families of soft recollisions 
come up naturally within the analytical i?-matrix theory 
of photoionization [19-24]. The ARM theory implements 
the intuition that the dominant effect of the ion’s electro¬ 
static potential on the photoelectron is a phase, and makes 
this intuition rigorous via the use of eikonal-Volkov wave- 
functions [34, 35] for the photoelectron. As in SFA treat¬ 
ments, this results in a time integral over the ionization 
time which, when approximated using saddle-point meth¬ 
ods, produces a quantum-orbit picture based on semiclas- 
sical trajectories. 

This first-principles approach has the advantage that it 
provides initial conditions for the trajectories by match¬ 
ing the eikonal-Volkov solution away from the core to the 
bound wavefunction near the core, using the WKB approx¬ 
imation for the bound wavefunction. In particular, the 
trajectories produced by this matching procedure are real¬ 
valued at the entrance to the classically-forbidden region, 
and this results in a nonzero imaginary position at real 
times, after exiting the classically-forbidden region. 

Thus, ARM operates with quantum orbits [36] which, 
in contrast to their purely real-valued counterparts in the 
Coulomb-corrected SFA [15-18], can change the ampli¬ 
tude of their contribution after ionization, i.e., after the 
tunnelling electron leaves the classically forbidden region 
[22,23]. Mathematically, this change in the amplitude 
comes through the imaginary part of the position, which 
results in an imaginary part of the action that directly af¬ 
fects the amplitude. Physically, this mirrors the dynamical 
focusing effect found in approaches that use the full classi¬ 
cal trajectory [25,26]. 

The importance and the effect of the imaginary part of 
the semiclassical trajectory is most strongly felt when the 
electron revisits the core. Mathematically, one needs to ex¬ 
tend the core potential into the complex plane. When the 
real part of the trajectory is small, this analytical continu¬ 
ation has a branch cut which cannot be integrated across. 
This branch cut must be managed carefully, as it can pre¬ 
clude several standard choices of integration contour, in¬ 
cluding in particular the contour going from the ioniza¬ 
tion time directly to the real time axis and then along 
it. In these cases, one needs a more flexible approach to¬ 
wards contour choice; to ease this choice we briefly describe 
user-friendly software [37, 38] to visualize the effects on the 
complex-valued position of different contour choices. 

Moreover, we present a method for consistently and pro¬ 
grammatically navigating the Coulomb branch cuts to cal¬ 
culate photoelectron spectra. This method is based on the 
fact that the branch cuts come in pairs which always con¬ 
tain between them a saddle point of the semiclassical dis¬ 
tance to the origin, \J r(f) 2 ; this saddle point is a solution 
to the closest-approach equation r (t) ■ v(t) = 0 on the com¬ 


plex plane. These times of closest approach offer a rich 
geometry of their own, both with complex quantum orbits 
and within the more restricted simple man’s model. More 
practically, by choosing appropriate sequences of closest- 
approach times, one can systematically choose correct in¬ 
tegration contours. 

Within this framework, soft recollisions appear as com¬ 
plex interactions between different sets of branch cuts, 
marked by close approaches between two or three closest- 
approach saddle points and by topological changes in the 
branch cut configuration. This happens only at very low 
transverse momenta, and is managed well by our method. 
Nevertheless, the increased time spent by the electron near 
the ion - mirrored in the ARM formalism by close ap¬ 
proaches between singularities and saddle points - results 
in an increase in the amplitude that reflects the photoelec¬ 
tron peaks seen in experiments. 

In the current approach, the Coulomb field of the ion is 
allowed to influence the phase of the wavefunction (and 
from there the ionization amplitude), but not the un¬ 
derlying trajectory. A first-principles approach based on 
the full trajectory is still lacking, but such a theory will 
likely require the trajectory to be real before it enters the 
classically-forbidden region during the ionization step. It 
is then likely that such a theory will contain many of the 
elements we describe for the laser-driven trajectories, in¬ 
cluding the branch cut behaviour and its navigation. In 
that sense, our work is a roadmap for those difficulties and 
their resolution. 

This paper is structured as follows. In Section I we ex¬ 
plore soft recollisions for the classical trajectories of the 
simple man’s model, giving simple approximations for the 
momenta that produce them, and we explore their scaling. 
In Section II we derive a simple version of ARM theory, 
emphasizing the features that give rise to complex-valued 
positions for the quantum orbits. 

We then examine, in Section III, the ways in which 
complex-valued positions give rise to branch cuts and the 
cases in which the latter make more sophisticated contour 
choices necessary. In Section IV we study the geometry 
of the times of closest approach, both within the simple- 
man’s model and for the ARM quantum-orbit picture, and 
explain their use for navigating integration contours around 
the Coulomb branch cuts. We then present the resulting 
photoelectron spectra in Section V and summarize our re¬ 
sults in Section VI. In addition, we provide supplementary 
information [39] with interactive 3D versions of several fig¬ 
ures in this paper. 


I. SOFT RECOLLISIONS IN CLASSICAL 
TRAJECTORIES 

The simplest approach to the dynamics of a photoelec¬ 
tron after tunnel ionization is known as the simple man’s 
model. In its quantum version it accounts for the tun¬ 
nelling process using the strong-field approximation, and 
then allows the electron to follow a laser-driven trajectory 
after it exits the tunnelling barrier. Within the SFA, the 
electron has velocity 


v(i) = P + A(t) 


(1) 



3 


where p is the canonical momentum registered at the de¬ 
tector. The electron is ‘born’ into the continuum at the 
origin at a complex starting time t s = f 0 + *tt which obeys 

2 (P + A(t s )) 2 + I P = 0, ( 2 ) 


where A(f) = —/F(t)di, F(t) = Fcos(wt)z is the laser’s 
electric field, which we take to be a monochromatic, lin¬ 
early polarized pulse, and I p = |ac 2 is the ionization po¬ 
tential. (We use atomic units unless otherwise noted.) The 
resulting trajectory is then 


r c i(t) = J (p + A(r))dr, (3) 

which is in general complex-valued. 

Throughout this work, we understand a ‘soft recollision’ 
to mean a trajectory which has a laser-driven turning point, 
with essentially zero velocity, in the neighbourhood of the 
ion [25, 26]. This can occur for a range of impact parame¬ 
ters, and the most relevant trajectories will pass within a 
few bohr of the ion. However, it is often easier to identify 
those laser-driven trajectories that ‘go through’ the ionic 
core, so we focus on them in the understanding that the 
neighbouring trajectories are the most relevant. 

Thus, to look for soft recollisions in the corresponding 
classical trajectories, we require that the real parts of both 
the velocity and the position vanish: 


Re [.Sci(tr)] = Re 


(j>z + A( t )) dr 


\Jt 3 


= 0 


v z (t r ) = p z + A(t r ) = 0. 


(4a) 

(4b) 


This can be expressed as 


F 

Zexit +Pz(t r - to) H — 2 (cos (wt r ) - cos(wf 0 )) = 0 (5a) 


UJ* 


F 

p z -sin (wt r ) = 0, (5b) 

uj 


where 


~exit — Re 


F 


J (p + A(r))dr 


= — cos(wto) (1 — COSh(tJTT)) 




( 6 ) 


is known as the tunnel exit. 

This system can be solved numerically, but it is more 
instructive to consider its linearized version with respect to 
p z , since all the soft recollisions happen at small energies 
with respect to U p . To do this we express the starting time 
as 
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FIG. 2. Scaling of the normalized momenta ujp z /F as a func¬ 
tion of the Keldysh parameter 7 for the first six soft-recollision 
trajectories. Dots show the exact solutions of (4) and lines show 
the linearized result ( 11 ). 


In the tunnelling limit of 7 <C 1, z e xit reduces to —I p /F as 
expected. 

The linearized system now reads 

f Pztr + ~2 (cos (ujt r ) - \J1 + 7 2 ) = 0 (9a) 

F 

p z -sin (ujt r ) = 0, (9b) 

w 

and to obtain a solution we must linearize t r with respect 
to p z . Examination of the numerical solutions of the orig¬ 
inal system shows that soft recollisions exist for every half 
period after uot r = 2 tt, so we write 


ujt r = (n + l)ir + u 5t r (10) 


with n = 1, 2, 3 ... indexing the solutions. With this we 
obtain from (9b) that 6t r ~ (—1 ) n+1 p z /F and cos (ojt r ) ~ 
(—l) n+1 . This gives in turn the drift momentum of the 
successive soft-recolliding trajectories as 


F v / 1 + 7 2 + (~ 1 ) 

W (n + 1)7T 


( 11 ) 


These are shown in Fig. 2, and are generally a good ap¬ 
proximation to the exact solution of the system (4). 

It is clear here that trajectories with odd and even n, 
which approach the origin from different sides, will behave 
very differently, and in particular will have very different 
scaling for low 7 . In particular, 

• for odd-n trajectories, in the tunnelling limit 


F q 2 7 K 

uj (n + l)27r (n + 1)27 t ’ 


( 12 ) 


to + * r T — t a 


— arcsin ( ^ (p z 
uj VF 



Pz 1 

F \Jl + 7 2 


— arcsinh ( 7 ), 

UJ 


where 7 = ujk/F is the Keldysh parameter, so that 


( 7 ) 


and the kinetic energy scales as y 2 for a fixed target 
species, whereas 

• for even-n trajectories 


P 


ST 

Z 


F 2 

UJ (n + l)7T ' 


(13) 


Zexit ~ (\/l +7 2 - l) ■ 


and the kinetic energy is a constant fraction of U p , as 
( 8 ) found previously [25,26]. This scales as 7 -2 . 
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These differences are grounded in much simpler scaling 
considerations. For even-?r trajectories, the centroid of the 
oscillation must advance past the ion by twice the quiver 
radius £ qu i v = F/uj 2 plus the tunnel length z ex it = \ r ) 2 z qu iv 
within the given number of laser periods, which implies a 
scaling of the form 

sr *^ / - q niv T 2-exit ^ 2 Ff UJ F . ^ , . 

P Z rri O / ^ V / 

1 Zir/LO uj 

and an energy which scales as 7 -2 . (Further, this sug¬ 
gests that the small observed deviations from this exponent 
[1,13] could also be explained by a shifted power law to ac¬ 
count for the tunnel exit, as confirmed in Ref. 40.) For the 
odd-ro trajectories, on the other hand, the centroid must 
only advance by the tunnel length in a multiple of a laser 
periods, and the scaling reflects this: 


-exit. _ Ip / F IpUJ 

~T~ ~ 2tt/w ~ ~F' 


(15) 


Thus, the fixed distance to cover in an increasing amount of 
time gives these trajectories their very low momenta at long 
wavelengths. Indeed, for argon in a 3.1 pm field of intensity 
10 14 W/cm 2 , as in Ref. 32, the highest odd-n momentum 
is p z « 0.02 a.u. This corresponds to a minimum energy of 
8 meV, which is consistent with the observed momentum 
spread of the ZES. 

The difference in scaling behaviour also suggests avenues 
of experimental research for testing this connection, as well 
as for probing the internal structure of the ZES. In partic¬ 
ular, it’s desirable to have experiments which increase the 
energy scale of the ZES and lift it above the consistent- 
with-zero level of the current experimental resolution, so 
that its details can be better examined. If it is caused by 
soft-recolliding trajectories of this nature, then its energy 
should scale as 7 2 I P , which points the way to experiment. 

This scaling is fairly unfavourable, since it is also neces¬ 
sary to keep 7 small to remain in the tunnelling limit. This 
means, then, that the ZES is best probed using high-/ p tar¬ 
gets, of which the most natural is Hc + ions - either as an 
ionic beam or prepared locally via sequential ionization or 
a separate pre-ionizing pulse. This brings a four-fold in¬ 
crease in I p compared to hydrogen, and allows for a wider 
range of intensities and wavelengths which keep 7 = kuj/F 
small but p s z r ~ 7K = k 2 uj/F large enough to be directly 
measured. 

Finally, we note that the momentum ratios between the 
two families of trajectories form universal sequences inde¬ 
pendent of the system’s parameters: 3/5, 5/7, 7/9,... for 
even -n trajectories and 1/2,2/3, 3/4,... for odd-n ones. 
The latter sequence should be difficult to resolve using cur¬ 
rent experimental setups but, if observed, would allow a 
richer window with which to observe the tunnel exit dis¬ 
tance [40]. 

It is also important to point out that, irrespective of 
the precise mechanism which translates the soft-recolliding 
trajectories into peaks in the photoelectron spectrum, this 
mechanism should apply equally well to both series of tra¬ 
jectories. As pointed out above, quantum-orbit approaches 
increase their associated amplitude via imaginary parts 
of the action, whereas classical-trajectory approaches give 
photoelectron peaks through dynamical focusing near the 


soft recollisions. In either case, the dynamical similarity 
between the two series of recollisions, evident in Fig. 1, 
should cause similar results for both. The reason the odd- 
n trajectories have not appeared explicitly so far seems to 
be only their very low energy . 1 


II. ANALYTICAL A-MATRIX THEORY 

We will now give a brief recount of ARM theory, as de¬ 
veloped in Refs. 19, 20, emphasizing the features that pro¬ 
vide initial conditions for the complex-valued quantum or¬ 
bit trajectories, thereby leading to branch cuts in the com¬ 
plex time plane. 

The i?-matrix approach was adapted from nuclear phy¬ 
sics [41] to model collisions of electrons with atoms [42] and 
molecules [43] , and it is designed to deal with electron cor¬ 
relation effects by confining full many-body dynamics to 
the interior of a sphere around the system, where they are 
most important, and using a single-active-electron approx¬ 
imation outside it. In a strong-field context, this permits a 
fuller account of multi-electron effects [44,45], and it also 
allows a smoother matching between the WKB asymptotics 
of the bound states near the ion with the perturbing effect 
of the ion’s potential on the continuum wavefunction on 
the outside region where this effect is small. 

The effect of the ion can be accounted for rigorously, in 
the outer R-matrix region, by using the eikonal approxi¬ 
mation [34,35] on the usual Volkov states, where the phase 
and amplitude of the wavefunction are approximated by 
the first terms in a semiclassical series in powers of H. This 
approximation yields continuum time-dependent wavefunc- 
tions of the form 

(r|p EVA (f)) = 1 e »(P+A(t))-r-4 //(p+A(r)) 2 dr 

(27t) 3 / 2 

X e —i//C/( r L(r;r.p,t))dr ) (jg) 

where U is the ion’s electrostatic potential, and 

r L (T;r,p,t) = r+ y (p + A(r')) dr' (17) 

is the laser-driven trajectory that starts at position r at 
time t and has canonical momentum p. These eikonal- 
Volkov states are (approximate) solutions of the time- 
dependent Schrodinger equation with the hamiltonian 

H = |p 2 + U[v) + r • F, (18) 

or, in other words, they are propagated in time as 

|p EVA (i)> = U(t,t') |p EVA (f')> for t > t' (19) 

where the propagator U(t,t') obeys the (approximate) 
Schrodinger equation 

d 

i—U(t,t ) = with U(t,t) = 1. (20) 


1 In addition, if the position of the tunnel exit is dropped from the 
analysis then all odd-n trajectories will degenerate to zero energy, 
which helps confound the situation. 
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The inner and outer regions are separated by a sphere 
around the ion of radius a. In strong fields it will be re¬ 
quired to satisfy the conditions 1 /k <C a <C I p /F , so it is 
well away from both ends of the tunnelling barrier; the 
wavefunction is correspondingly split into an inner and 
outer wavefunctions. This carries the problem that the 
hamiltonian is no longer hermitian in either region (as the 
usual integration-by-parts proof leaves uncancelled bound¬ 
ary terms), so the resulting system has two non-hermitian 
Schrodinger equations coupled through their boundary con¬ 
ditions. 

The problem is rigorously addressed by introducing a 
‘hermitian completion’ of the hamiltonian, by means of the 
Bloch operator 


L = S(r — a) (’ (21) 

for which H + L is hermitian in the inner region and H — L 
is hermitian in the outer region, with the Bloch operator 
cancelling out the boundary terms in the usual integration 
by parts. The system can then be re-expressed as two cou¬ 
pled, hermitian, inhomogeneous Schrodinger equations: 

t) = (H + T)4' in (r, t) - L^ 0 ut (r, t) (22a) 

*^^out(r, t) = (H - L)^ out { r, t ) + TT in (r, t). (22b) 

at 

In particular, since the Bloch term LT is local to the 
boundary, we can implement the boundary conditions by 
substituting the wavefunction from the other side of the 
boundary in both equations, which then acts as a source 
term. The solutions to the coupled equations then auto¬ 
matically obey the boundary conditions. 

Once with a well-formulated system, one can apply the 
relevant approximations. In practice, the equation for the 
inner region is not affected as long as the ionization rate is 
not too great, since then the Bloch term L\P out can be ne¬ 
glected, or its main effects can be modelled by incorporat¬ 
ing a ground-state depletion factor; furthermore, by an ap¬ 
propriate choice of the constant b in the Bloch operator, the 
ground state can be made an eigenstate of L. (In partic¬ 
ular, choosing b= 1/n implies that Ld/ g = —n5{r — a)\& fl .) 
Similarly, the Bloch term in the outer region’s hamiltonian 
can be neglected as the ionized wavefunction will generally 
be far from the boundary. 

Given these approximations, one can then write down the 
formal solution to the outer Schrodinger equation (22b) as 

1 9(t)) = -i[ U(t,t')L\% n (t'))dt'. (23) 

J — OO 


The quantity of interest is the photoelectron momentum 
amplitude at a time T long after the pulse is finished, which 
is then 


o(p) 


(p|tf(T)) = -i [ T (p\U(T,t)L\* in (t))dt 

J — OO 

-i [ T (p EVA (t)\L\* g (t))dt. (24) 

J — OO 


To calculate the spatial matrix element, we work in the 
position basis, so 


*(P) = ~i fdv fdt (p EVA (f)|r> <r|L|* fl (i)> 

J J —OO 

= —i f clr f dt (-K.)S(r — a) (p + A(t)|r) (r|^ s ) 

J J —OO 

x e iI P t +1 /y(p+A(r)) 2 dr e i/y [/(r L (r;r,p,t))dr 

=—i f ^e i/pt+ 5 /tIp+A-C)) 2 ^ 

J —OO 

x (— k) J dr< 5(r — a) (p + A(t)|r) (r|^ 3 ) 
x e */r c/ ( rL ( r ; r ’P’ t )) dT . (25) 


This calculation produces, then, an ionization amplitude 
a(p) in the form of a temporal integration of a phase factor 


giSV(i) = exp 


° t+ \J T (P + A(r)) 2 dr 


(26) 


a Coulomb phase factor arising from the integration of the 
ionic potential along the laser-driven trajectory, and a spa¬ 
tial factor involving a Fourier transform of the ground state. 
This amplitude is structurally similar to the SFA ampli¬ 
tude: in particular, the phase factor (26) is present in the 
SFA. It oscillates at the frequencies I p and U p , which are 
fast compared to the laser-cycle timescales that the inte¬ 
gration takes place on, and this makes the integral highly 
oscillatory. This enables, in turn, the application of the 
saddle-point approximation [46]. However, this is compli¬ 
cated by the fact that the eikonal Coulomb correction to 
the Volkov phase, 


exp 


i J U (r L (r; r, p, t))dr 


(27) 


couples the spatial and temporal integrations in a compli¬ 
cated fashion. This factor, moreover, is crucial in ensuring 
that the resulting amplitude is independent of the (unphys¬ 
ical) boundary radius a. 

To disentangle both integration steps, one first performs 
the temporal saddle-point integration, with an r-dependent 
result since the phase now includes a spatial term e lA (*b r . 
The resultant saddle point t a is generally not far from the 
standard SFA saddle t s = t 0 + itt, which obeys 

p (P + A(f s )) 2 + I p = 0. (28) 

Moreover, the variations in t = t a in the Coulomb correc¬ 
tion term (27) exactly cancel out the variation with respect 
to a in the ground-state wavefunction [see 19, 20, 24]. 

The result of this matching procedure [19, 20, 24] is a 
simple expression for the ionization amplitude with most 
factors evaluated at the complex saddle-point time t s of 
(28), 


o( p) =e i/ ^e“ 2 ^ (p+A(T))2dT 

x e ~iJu U (Z P+ A ^') dT ') dT i?( p ) ; (29) 


where t K = t s — i/n 2 and the final amplitude will be a 
sum over all saddle points t s with positive imaginary part. 
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(Hereafter, however, we focus on the contribution from a 
single half-cycle.) In this expression the spatial factors have 
been incorporated into a shape factor 


R( P) 


ina 2 

y/iSy(t a ) 


— p~ i ^ta U (ft, P+A(r')dr')dr 

2tt 6 


x e~^ +A ^ r (r|M/ g ) | r=a , (30) 


which is independent of a [19,20,24], This shape factor 
is essentially the Fourier transform of the ground state 
over the spherical boundary; in this sense it is analogous 
to the planar boundary used in partial Fourier transform 
approaches [47]. 

The ARM ionization amplitude Eq. (29) is similar to the 
SFA result, and it admits a clear quantum-orbit trajectory 
interpretation [36]. The electron is ionized at time t s and 
propagates to its detection, at real and large time T, with 
velocity p + A(t) and position 


r c i (t) = 


(p + A(r)) dr. 


(31) 


correction factor 


exp 


-if U (r c i(r))dr , 
Jt n 


(32) 


from a time t K = t s —i/n 2 just under the saddle point. This 
shift arises here rigorously as a result of matching the outer- 
region wavefunction to the WKB asymptote of the bound 
wavefunction. Since its origin is intrinsically spatial, it is 
preferable to other regularization procedures such as the 
ionization-rate-based scheme of Ref. 13. 

Finally, it is important to note that the process that 
leads to the Coulomb correction factor (32) is little more 
than the analytic continuation of the original factor (27), 
which appeared in the initial real-valued temporal integral 
(25). This uniquely determines the trajectory to use for 
the Coulomb correction, with the initial conditions set by 
the matching procedure. The trajectory r c i (t), along which 
the Coulomb correction is integrated, obeys the laser-only 
equation of motion, 


fd(t) = —F(t), (33) 


Along the way, it acquires ‘phases’ given by the exponenti¬ 
ated action of the different relevant energies: e zIpts from the 
bound energy before ionization, e _5 /t s (p+ A ( T b dr f r0 m the 


kinetic energy after ionization, and c/ ( rcl O) dT f rom 

the Coulomb interaction. It is important to remark that 
these factors are no longer pure phases since the exponents 
are now complex, which greatly affects their amplitude and 
indeed completely encodes the (un)likelihood of the tun¬ 
nelling process. 

The complex integrals in Eq. (29) are complex contour 
integrals which can be taken, in principle, over any con¬ 
tour which joins both endpoints. In practice, the stan¬ 
dard contour goes directly down from t s to its real part 
to and then along the real axis until the detection time T, 
as shown in Fig. 3. This conveniently separates the imagi¬ 
nary part of the contour, which yields an amplitude reduc- 

that encodes the tunnelling 


tion 


e Upta g-1 / t ° (P+A(r)) 2 dr 


probability, and the normal propagation along the real time 
axis. As we shall see later, when the Coulomb interaction 
is included over one or more laser periods, this choice of 
the second leg of the contour is no longer necessarily ideal 
or even allowed. 

The main effect of the matching procedure is that the 
r-dependent laser-driven trajectory (17) has been replaced 
by a classical trajectory of Eq. (31) starts at the origin, 
but this trajectory is only integrated over, in the Coulomb 
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FIG. 3. The time contour for integrating the exponents in 
Eq. (29) is usually taken from t s down to its real part to and 
then along the real axis until the detection time T. 


as opposed to the full equation (including both the laser 
and the Coulomb field) used in the CCSFA approach [15- 
18]. This is a direct consequence of the appearance of r c i 
in the EVA wavefunction (16). Work to include the effect 
of the ion on the trajectory should ideally be directed at 
its inclusion at the level of this wavefunction. 

The initial conditions are also different to those used in 
the CCSFA approach. In ARM theory, the initial position 
is obtained from first principles and it is real at the entrance 
of the tunnelling barrier instead of its exit, as is imposed 
externally on the CCSFA trajectory. The consequences of 
this can be dramatic, as we show below, and among other 
things they preclude certain choices of integration contour 
in (32). In this light, calculating semiclassical trajectories 
for the full equation of motion with the correct initial con¬ 
ditions in complex time and space appears to be a far more 
difficult task. 

In the next section we describe the origin of these dif¬ 
ficulties, the mathematical challenges that arise, and how 
they can be addressed within the ARM theory. 


III. EMERGENCE OF TEMPORAL BRANCH 
CUTS 

The analytical i?-nratrix theory allows us, then, a tra¬ 
jectory-based approach for accounting for the effect of the 
ionic potential on the photoelectron. Here the potential 
is integrated over the electron’s complex-valued trajectory 
(31) to give the correction factor (27), which now has a non¬ 
trivial amplitude. The trajectory starts at a real position, 
r c i(t K ) ~ — z/re, but from then onwards both the integrand 
and the integration time are complex, so the position is 
therefore complex. 

The imaginary part of the position has two main sources. 
The first is along the transverse direction, since any amount 
of real-valued transverse momentum will accumulate an 
imaginary transverse position ipj_Tr between the complex¬ 
valued t s = to + itt and any real time t, from the ‘tunnel 
exit’ t = to onwards. The longitudinal position, on the 
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FIG. 4. Semiclassic.al trajectories on the complex 2 plane, where 
t covers the ‘under-the-barrier’ segment from t K to to, shown in 
the inset. When p z = 0 the trajectory is completely real, but for 
nonzero momenta the imaginary part is increasingly important. 
Paths are labelled by the normalized longitudinal momentum 
top z /F. Hereafter k = 1.07 and F = 0.05, corresponding to 
argon in a 9 x 10 13 W/cm 2 field, and A = 1pm (so 7 = 0.98) 
unless otherwise noted. 


other hand, also obtains imaginary components through 
the integral of the vector potential, which is now explicitly 
complex valued, so the end result is the nonlinear depen¬ 
dence shown in Fig. 4. 

Since the integrand in (31) is real for real times, this 
imaginary position will not change along the real time axis, 
so it will remain present until the electron is detected at a 
large time T, at which the real position will be much larger. 
During a recollision, however, the real position becomes 
small. At this time, the imaginary coordinate dominates 
the position, which becomes an issue when calculating the 
ion’s electrostatic potential. The simplest choice for that 
is the Coulomb potential, 


U( r) 


1 _ 1 
Vr 2 \/ x 2 + y 2 + z 2 ’ 


(34) 


which is analytical in this form. Other potentials are also 
tractable, as long as they are analytical functions on the 
real axis, but pose similar mathematical challenges. 2 For 
the Coulomb potential, the extension procedure requires 


2 In particular, any atomic or molecular charge density is accept¬ 
able as long as an explicit analytical formula is available, which in 
practice is excessively restrictive. Extended spherically-symmetric 
charge densities like exponentials or softened Coulomb potentials 
exhibit similar behaviour but with the branch cut pushed back 
by the characteristic width of the distribution. In contrast, the 
potential created by a gaussian charge shows drastically different 
behaviour, with the branch cuts replaced by growth ase +|r| where 
Im(r 2 ) <0. If the ionic charge distribution is only known numer¬ 
ically - e.g. via quantum chemical methods - then both models 


the use of a square root, which in general poses at most 
an integrable singularity, but for complex arguments it has 
a branch cut along the ray r 2 £ (— 00 , 0) which cannot be 
ignored. 

This branch cut is reached, precisely, on recollisions 
where the imaginary part dominates a nonzero real part. 
Indeed, the squared position can be decomposed as 

r 2 = Re(r) 2 — Im(r) 2 + 2 i Re(r) • Im(r), (35) 

which means that, when Im(r) 2 > Re(r) 2 , the trajectory 
will reach the branch cut at the point where the projection 
Re(r) • Im(r) changes sign. This is marked by a discontin¬ 
uous change in the square root Vr 2 , with a sudden change 
in the sign of the imaginary part Im ^ x/r 2 ^. If one inte¬ 
grates across this discontinuity, the dependence on t of the 
Coulomb integral (27) ceases to be analytic, and one loses 
the freedom that allowed the real contour of (25) to be 
deformed to pass through the saddle point t s in the first 
place, and the results are no longer meaningful. 

The solution to this problem is to view the entire poten¬ 
tial as a function of time, 

UMt)) =~viw (36) 

as a single analytical function of the complex-valued vari¬ 
able t. The branch cuts in U are imprinted on the time 
plane via the conformal mapping t 1-7 r c i(t) 2 . We show in 
Fig. 5(a) an example of the Coulomb potential’s behaviour, 
as a contour map of the function \/r c \(t) 2 (which has the 
branch cut structure of U but omits its singularities). The 
essential features of this function are the branch cuts, which 
are sketched in Fig. 5(b) in red. In particular, the standard 
integration contour of Fig. 3 - straight down from t s and 
then along the real axis, shown dotted in Fig. 5(b) - can 
indeed cross the branch cuts when the electron returns near 
the ion. 

This means that to preserve the analyticity of (27) one 
must deform the integration contour away from the real 
time axis until the integrand is continuous and analytic. 
This will correspondingly change the way the complex po¬ 
sition r c i(i) moves in space, with the main effect of min¬ 
imizing the imaginary part of the position at the time of 
recollision; this is exemplified in Fig. 6. 

The relationship between the chosen temporal con¬ 
tour and the corresponding trajectory in complex position 
space, particularly along z, is in general complicated and 
hard to visualize. Similarly, the final momentum p has a 
strong effect on the potential U(r c \(t)) as a function of time, 
sometimes very sensitively (such as near a soft recollision). 
These variables are all intertwined, with the momentum de¬ 
termining the branch cut structure and therefore the pos¬ 
sible contours on both the time and space complex planes. 

To help disentangle these relationships, we have devel¬ 
oped a software package [37] which enables real-time visu¬ 
alization of the effect of a particular choice of contour and 


pose conflicting demands, and numerical evaluation via integration 
of a Coulomb kernel fails to give analytical results. We will analyse 
these differences in detail in a future publication. 
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FIG. 5. (a) Contour plot of the complex distance to the origin, ^/r c i(i) 2 , with the coloured background along the real part and 

the red, black and blue orthogonal lines representing, respectively, positive, zero and negative imaginary parts. The white lines are 
branch cuts where the real part is zero and the imaginary part discontinuously changes sign. In (b) we show a simplification of 
this picture, with the red lines showing the branch cuts, and the thin green lines the positive real axis of r 2 . The shaded regions 
indicate the complex times for which Re(r c i(t) 2 ) is negative, which are undesirable when using gaussian and numerically-obtained 
ionic potentials, as they would be unphysical there. The momentum displayed, p = (0.02,0, 0.8), is such that the standard contour 
(integrating down to the real part to of the saddle point t a and then along the real axis, shown dotted in (b)) crosses a branch 
cut. Instead, one should choose a contour which passes through a time between the branch cuts, which we label I C a and explore in 
detail in Sec. IV. 


momentum on the space trajectory and on the temporal 
branch cut structure. This software is described in detail 
in Ref. 38 and we encourage the reader to use its visualiza¬ 
tion tools to explore the implications of our results. 


IV. TIMES OF CLOSEST APPROACH 

We see, then, that the existence of branch cuts which 
can cross the real axis can preclude the use of the standard 
integration contour for the integral in (27), and one needs 
to choose a contour which passes through the ‘slalom gate’ 
left by the branch cuts of U{r c \{t)) as in Fig. 5(a). If only 
a few momenta are involved then this can be done on a 
case-by-case basis, but the computation of photoelectron 
spectra requires a more programmatic approach, which is 
able to automatically choose the correct contour for each 
given momentum. 

To obtain this approach, we examine in detail the space 
between the two branch cuts, as shown in Fig. 7 . Each 
branch cut is a contour of constant Re(y / ryi(t) 2 ) = 0, 
which means that the neighbouring contours must closely 
follow its direction, and circle around it when it terminates 
at the branch point. In the ‘slalom gate’ configuration of 


Fig. 5(a), the only way to combine the curvatures of the 
contours near the facing branch cuts is to have a saddle 
point in between them. This saddle point is the crucial ob¬ 
ject which enables the automatic choice of a contour which 
avoids the branch cuts, as it is ideally situated between 
them. Moreover, it has a deep physical and geometrical 
significance, which we explore in this section. 

It is clear in Fig. 7 that any contour which crosses the 
complex t plane from left to right must have a real distance 
to the origin Re(-y/r c i(f) 2 ) which decreases, reaches a min¬ 
imum, and then increases again. For contours which cross 
the branch cut, this minimum is zero, and it is accompa¬ 
nied by a discontinuous change in the imaginary part. For a 
path crossing through the saddle point, on the other hand, 
the minimum of Re(-\/rci(i) 2 ) is at its maximum value, and 
this minimum occurs precisely at the saddle point. 

For this reason, we call the saddle point the time of clos¬ 
est approach, and we label it t CA . To be precise, then, a 
contour that passes through t = t CA maximizes the mini¬ 
mum value of the real part of the distance to the origin, 
Re(\/r c i(f) 2 ). Intuitively, it permits the furthest possible 
approach to the ion. 

Similarly, when calculating the Coulomb potential along 
such a contour, this choice minimizes the maximum value 
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Re(z) 


FIG. 6. Trajectories in the complex 2 : plane corresponding to 
the contours shown in Fig. 5(b). The trajectory starts at 2 = 
— 1/k at time t K , and departs somewhat from the real position 
axis as the time goes down to the real axis at to, as shown in 
Fig. 4. Along the real axis, on the standard contour shown 
dashed, the electron goes to large negative Re( 2 ) before turning 
around towards the core. It then reaches the ion with a large 
imaginary part, which causes a discontinuous jump in the square 
root of (34). Deforming the contours to avoid the branch cuts, 
shown in the solid line, minimizes the imaginary part of 2 at the 
moment of recollision; it is then slowly regained before detection 
at a large real time T. The closest approach time t C A marks the 
minimum value of Re(r c i(f) 2 ) once the x coordinate is taken 
into account; there 2 c i is small but nonzero. 


tines. 

To find the saddle points, one simply looks for zeros of 
the derivative of i/r c i(t) 2 . This can be further reduced to 
the zeros of gg [r c i(f) 2 ], so the criterion is simply 

M*ca) • v(f CA ) = 0. (37) 

This equation is deceptively simple, and one must remem¬ 
ber that the left-hand side is a complex-valued function of 
time through Eq. (31). Nevertheless, it has a compelling 
physical interpretation, for if a classical electron passes near 
the nucleus then it is closest to the origin when its velocity 
and its position vector are orthogonal. 

In this spirit, then, it is worthwhile to investigate the 
classical solutions of (37) before exploring the solutions in 
the complex quantum domain. As we shall see, both do¬ 
mains exhibit rich geometrical structures which are closely 
related to each other and, moreover, the times and mo¬ 
menta of soft recollisions emerge naturally as crucial ge¬ 
ometrical points within both structures. After exploring 
the geometrical implications in both contexts, we shall use 
this knowledge to automatically generate correct integra¬ 
tion contours for any momentum. 


A. Classical times of closest approach 


of the real part of 1 />/r c i(t) 2 and its absolute value (for 
valid contours which do not cross the branch cuts), so that 
the Coulomb interaction is kept as bounded as possible. As 
long as the potential [/(r cl (f)) stays continuous and ana¬ 
lytic, this is not essential, as the integral in (27) does not 
change. However, this contour optimizes the applicability 
of the approximations that led to (27), and it admits the 
clearest physical interpretation by keeping the imaginary 
part of the position within the tightest bound possible at 
the points where this is necessary. Additionally, it min¬ 
imizes the difficulties faced by numerical integration rou- 
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FIG. 7. Closer view of the region between the branch cuts 
in Fig. 5(a). The space between any two branch cuts always 
contains a saddle point t C A, as this is the only way to join the 
curvatures of the contours next to each branch cut. For contours 
that pass through it, the saddle point marks the minimum of 
Re ^y/r c i (t) 2 J . For valid contours that do not pass through t c a, 
this minimum will be smaller. 


The classical closest-approach times are naturally 
defined via 

Re [r cl (ip") • v(t c c ‘“)] = Re MO] ’ v (*cD = 0. (38) 

The problem in using them as a tool, however, is that both 
of their defining equations, (37) and (38), have several other 
solutions, and these are not always distinguishable from the 
desired closest-approach times. In particular, the turning 
points of the classical trajectory, away from the core, are 
also solutions of (37), since they are also extrema of r 2 . 
Thus, for instance, the surface in Fig. 5(a) contains saddle 
points at uit ~ 7 t/4 and cot ~ 37 t/4 , which correspond to 
the turning points shown in Fig. 6 to the right and left of 
the position at t 0 , respectively. Thus, to be able to use 
the closest-approach times as an effective tool to avoid the 
branch cuts one needs to distinguish the crucial mid-gate 
t CA points from the other solutions. 

In certain cases this is easy, such as for on-axis collisions 
with p± = 0. Here (38) reads 

Re [Mc'D] McD = °> ( 39 ) 

and it solutions cleanly separate into turning points, with 
v z{tc“) = 0, and collisions with z c i(i) = 0, as shown in 
Fig. 8(a). The turning points can additionally be classified 
as minima and maxima of r c \(t) 2 , shown respectively in 
green and red, by evaluating the sign of ^2-r c i(t) 2 . At 
nonzero p±, it is the collisions that will turn into useful 
closest-approach times. 

Within this context, the classical soft recollisions stud¬ 
ied in Section I appear naturally as the intersection be¬ 
tween these two curves: the times when both v z (t) = 0 
and z c \(t) = 0. Moreover, at these intersections the num¬ 
ber of available roots changes. Thus, as p z is swept up 
across the intersection marked a in Fig. 8(a), an inward 
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turning point turns into an outward turning point flanked 
by two closest-approach points. The classical trajectory 
exhibits exactly this behaviour, as shown in Fig. 9(a), with 
this change in p z . 

The roots of (39) can also merge at the extremes of 
the sinusoidal turning-point curve of Fig. 8(a). At these 
points, the longitudinal momentum becomes greater than 
the oscillation amplitude F/w, and the velocity v z (t) = 
p z — sin(wt) no longer changes sign. The resulting be¬ 
haviour of the trajectories is shown in Figs. 9(b) and (c), 
and resembles pulling a winding string until the turns are 
straight. 

The closest-approach solutions of (38) become more in¬ 
teresting when one allows a nonzero transverse momen¬ 
tum p x . Here the solutions form a single coherent surface, 
shown in Fig. 8(b), that consists of a number of bounded 
lobes joined together at the soft recollisions, which locally 
look like cones. Thus, it is possible to continuously con¬ 
nect any two roots of (38) via a path on the surface. More 
specifically, the inward turning points and the recollisions 
shown in green in Fig. 8(a) can be smoothly connected via 
the p x ^ 0 component of the surface, which precludes the 
existence of a simple criterion to distinguish one from the 
other in the general case. 

On the other hand, the outward turning points can still 
be distinguished, as they are local maxima of r c i(t) 2 . These 



0 7T/2 n 3 ttI 2 2n 5 ttI2 3jt 7ji/2 4n 


nl2 


3tt I 2 


5tt/2 


In 12 


(a) 15 

i. 

4. 0.5 

-0.5 

- 1 . 

0 


oA 


FIG. 8. The classical closest-approach times with p± = 0, 
satisfying (39), separate into two curves (a): turning points 
with v z = 0 and collisions with z c i = 0 , with their intersections 
representing soft recollisions. At the points marked a, b and c 
the different roots merge and disappear, with the corresponding 
trajectories displayed in Fig. 9. For nonzero p±, the solutions 
of the vector equation (38) form a single coherent surface (b) 
with a sequence of bounded lobes which connect at the soft 
recollisions, where the surface is locally a cone. Local minima 
and maxima of r c i(t) 2 are shown respectively in green and red in 
both panels. At the boundary between the two, a maximum and 
minimum meet, merge and disappear, as shown in the trajectory 
of Fig. 9(d); the point then leaves the surface. The side and 
bottom panels show the projections of the surface on p z and p x 
respectively. An interactive 3D version of this figure is available 
in the Supplementary Information [39]. 
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FIG. 9. Classical trajectories from near the critical points 
marked a-d in Fig. 8, where the closest-approach roots of (38) 
can merge and disappear, indexed by their reduced momentum 
( ujp x /F,iop z /F ). The dots show closest-approach points, which 
satisfy (38) and for which the tangent to the trajectory is or¬ 
thogonal to the radius vector. In the neighbourhood of a soft 
recollision, (a), an inward turning point turns into an outward 
turning point flanked by closest-approach points as the momen¬ 
tum increases. In (b) and (c) two turning points, a maximum 
and a minimum of r c i(f) 2 , merge and disappear as p z goes past 
the oscillation amplitude F/uj, either in the positive (b) or neg¬ 
ative (c) direction, so v z no longer changes sign and no turning 
points occur. Similarly, as the transverse momentum p x in¬ 
creases in (d) the x component of the velocity becomes too great 
for the tangent to be orthogonal to the radius vector; at that 
point a maximum and minimum of r c i(f) 2 merge and disappear, 
leaving r c i(£) 2 to grow monotonically. 


maxima are shown in red in Figs. 8(a) and (b), and they 
form the “left-facing” side of the surface in Fig. 8(b). 
More specifically, any horizontal line of constant momen¬ 
tum must enter the surface through a maximum (red) and 
leave it through a minimum (green), because the minima 
and maxima must alternate for any given trajectory. Thus, 
the red (maximum) side of the surface points towards neg¬ 
ative t. and the green (minimum) side points towards pos¬ 
itive t. 

At the boundary between both parts of the surface, a 
maximum and a minimum merge and disappear, and the 
trajectory will then behave as shown in Fig. 9(b), (c) or (d), 
depending on which direction the boundary is approached 
(i.e. towards positive p z , negative p z , and increasing p x , 
respectively). Horizontal lines of constant momentum will 
be tangent to the surface at this boundary, and the corre¬ 
sponding trajectory will have a double root of (38). 


B. The quantum solutions 

The quantum solutions have a richer geometry, with an 
additional dimension, imaginary time, to be occupied. The 
main effect of this is to increase the number of available 
solutions: while in the classical case two real solutions of 
(38) can merge and disappear, in the quantum case the 
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complex solutions of (37) are not lost, but will move into 
imaginary time and remain present. 

In general, the quantum solutions will be close to the 
classical ones when the latter exist. As one approaches the 
end of a lobe on the classical solution surface of Fig. 8(b), 
however, the quantum solutions approach each other close 
to the real axis and then diverge into positive and negative 
imaginary time, keeping a relatively constant real part. If 
one then projects this to real times, the result is a pair 
of surfaces which closely follow the red and green parts of 
the classical surface, and then diverge into roughly parallel 
planes as they reach the boundary. We show this behaviour 
in Fig. 10. 

The first few closest-approach solutions are relatively 
easy to handle, and depend smoothly on the momentum. 
This includes the first minimum and maximum of r c i(f) 2 , 
like the ones in Fig. 5(a), the birth time t s itself, and a con¬ 
jugate solution with negative imaginary part which should 
be ignored. These solutions occupy specific regions of the 
complex t plane, as shown in Fig. 11(a), and can be con¬ 
sistently identified. Moreover, these solutions exhibit close 
approaches at tot ss 7r/2 and tot ss 3n/2 which are the 
quantum counterparts of the classical maximum-minimum 
mergers shown in Figs. 9(b) and (c). These are evident in 
Fig. 10 as the converging surfaces at those times, and are 
of relatively limited interest. 

The most important t CA close approaches occur at and 
near the soft recollisions, shown inside the gray rectangle of 
Fig. 11(a) and in Fig. 11(c), with a complicated momentum 
dependence which we explore below. In the quantum do¬ 
main, soft recollisions again represent interactions between 
three different closest-approach roots. Unlike the classical 
domain, however, the roots do not merge; instead, two of 
them move into imaginary time after a three-way avoided 
collision, shown in Fig. 11(c) between the points marked 9 
and 10. The proximity between the multiple saddle points 
marks the increased time the electron spends near the ion 
in the neighbourhood of a soft recollision. 

More interestingly, this three-way collision marks a cru¬ 
cial topological change in the configuration of the branch 
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FIG. 10. The quantum solutions of the closest-approach equa¬ 
tion (37) form multiple surfaces which wrap around the clas¬ 
sical solutions of Fig. 8(b), closely following the lobes where 
they exist and departing at the edges to form pairs of parallel 
surfaces with imaginary parts of opposite sign. Black dots rep¬ 
resent largely real solutions, with red (blue) dots representing 
solutions with positive (negative) imaginary part. An interac¬ 
tive 3D version of this figure is available in the Supplementary 
Information [39]. 
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FIG. 11. The quantum closest-approach times, which satisfy 
(37), on the complex time plane. In (a) we show the closest 
approach times corresponding to a grid in momentum space 
(shown inset). These are generally grid-like in the time plane, 
though after tot = 3n/2 some solutions lie close to the real 
axis and cannot be discerned in this view. At the cusps of 
Figs. 8(b) and 10, which correspond to soft recollisions, the reg¬ 
ular grid-like behaviour breaks and the solutions can no longer 
be uniquely tagged. This can be seen by following the semicir¬ 
cular path in momentum space shown in (b), for which there 
are three solutions in the gray rectangle of (a), shown in de¬ 
tail in (c). Going once around the semicircle moves the fcA 
along the bow-shaped curves and, upon returning to the initial 
point, permutes them cyclically. In particular, this path con¬ 
tains an avoided collision between the points marked 2 and 3, 
and a three-way avoided collision between points 9 and 10. This 
three-way interaction marks the soft recollision itself. An inter¬ 
active 3D version of this figure, with considerable additional 
detail, is available in the Supplementary Information [39]. 


cuts associated with the recollision, as shown in Fig. 12. 
Each of the outer saddle points, tc A and t®, has a pair 
of branch cuts associated with it, in the same ‘slalom gate’ 
configuration as in Fig. 7, and these go off into imaginary 
time. However, the way in which they do so changes as 
the longitudinal momentum p z passes the momentum pf 
of the soft recollision. 

For p z below pf , as in Fig. 12(a), the branch cuts loop 
back to imaginary time without crossing the real time axis. 
As with the low-momentum trajectory of Fig. 9(a), the 
trajectory does not quite reach the collision, and the asso¬ 
ciated branch cuts do not force a change of contour. At 
p z = pf, however, the branch cuts touch and reconnect, 
and for p z > pf the topology changes to the one shown in 
Fig. 12(b). Here the trajectory does pass the core, and the 
associated branch cuts do cross the real axis, forcing the 
integration contour to change and pass through the gates. 
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This process has profound implications for the ionization 
amplitude, because these drastic changes in the integrand 
occur precisely when it is largest. Thus, choosing the wrong 
contour in this region accounts for the largest contributions 
to the integrand, with a correspondingly large effect on the 
integral. More surprisingly, once the contour is forced to 
pass through the ‘gate’ t CA s, for p z just above pf, their 
contributions have the effect of suppressing the ionization 
amplitude there. 

To see how this comes about, consider the integral 
fU(r c i(t))dt for the configuration of Fig. 12(a). Here 
y/i'd(t) 2 has a minimum at the central saddle point, tcH , 
and this translates into a maximum of l/\/ r c i(i) 2 which 
dominates the integral. At this point, the approach dis¬ 
tance r* = \J r c i(icA ) 2 is dominated by a modest and pos¬ 
itive imaginary part. This means that £/* = — 1 /r* is large 
and imaginary, and therefore the correction factor e ~ l I udt 
has a large amplitude. 

On the other hand, in the configuration of Fig. 12(b) 
the integral is dominated by the ‘gate’ closest-approach 

times, for which r* = \J r c i(fc* ) 2 ~ \J r c i(^cA ) 2 is mostly 
real and much smaller than r*. The corresponding po¬ 
tential t/' = — 1/r* is then large, real and negative, and 
—ill* is along +i. However, here the line element d t must 
slope upwards with a positive imaginary part to empha¬ 
size the contribution of the saddle point, and this then 
gives — i f U(r c \(t))dt a large and negative real part. This, 

in turn, suppresses the amplitude of the correction factor 
g —i f Udt 

This effect is visible in photoelectron spectra as a large 
peak just below the soft recollision, followed by a deep, 
narrow dip. In an experimental setting, the dip will almost 
certainly get washed out by nearby contributions unless 
specific steps are taken to prevent this, but the peak will 
remain. In addition, this effect mirrors the redistribution 
of population seen in classical-trajectory based approaches, 
where the peaks caused by dynamical focusing represent 
trajectories taken from other asymptotic momenta, whose 
amplitude is reduced. 

We now turn to the momentum dependence of the 
closest-approach times near the soft recollision, which again 
presents interesting topological features. The main prob¬ 
lem is illustrated in Figs. 11(b) and (c): the different solu¬ 
tions of (37) mix, and there is no longer any way to distin¬ 
guish them from each other, as there was in the classical 
case. More specifically, traversing a closed loop in momen¬ 
tum space, like the semicircle shown in Fig. 11(b), will move 
the roots around in such a way that when one returns to 
the initial point the overall configuration is the same, but 
the saddle points have been permuted cyclically. 

Topologically, this means that the surface defined by (37) 
(a two-dimensional surface in a four-dimensional space) 
does not separate into distinct components; instead, the 
surface has a single connected component after uit = 3ir/2. 
On the other hand, the surface itself remains singly con¬ 
nected. Both of these behaviours are explored in Figs. S3 
and S4 in the Supplementary Information [39]. 

This mixing behaviour is unusual in the quantum orbit 
formalism, where the norm is for rather elaborate indexing 
schemes to be possible [28,48], partly because there is usu¬ 
ally a single free parameter that governs the motion of the 



(c) v(t) 2 (d) v(t) 2 



6.1 6.2 6.3 6.4 6.5 6.6 6.1 6.2 6.3 6.4 6.5 6.6 


(e) Branch cut sketch (£) Branch cut sketch 



FIG. 12. During a soft recollision, the quantum times of clos¬ 
est approach will perform a three-way close approach, like that 
shown in Fig. 11(c) at the first soft recollision between the points 
9 and 10. At this close approach, the branch cuts associated 
with these saddle points will reconnect and change topologies, 
as shown in (a) and (b). This happens near, but not at, the 
classical soft-recollision momentum studied in Section I. At the 
point of the topological change, the outer saddlepoints and 

f o\ 

tcA emerge from the classically forbidden region, and the real 
part of their kinetic energy Re(^v(t) 2 ) changes sign, as shown 
in (c) and (d). After the change, the outer saddlepoints do not 
require tunnelling to get to, and should be crossed by the inte¬ 
gration contour. We show in (e) and (f) a sketch of the relevant 
branch cuts, analogous to Fig. 5(b), along with the appropri¬ 
ate integration contour for each topology. In (g) and (h) we 
show the corresponding trajectories in the complex z plane for 
those contours. The motion along 2 is similar to Fig. 9(a), with 
a modest imaginary part as in Fig. 4, and if the core is not 
reached does not represent a problem. For the slightly higher 
momentum of (h), however, the trajectory passes the origin and 
would cross the associated branch cut if taken along the dashed 
integration contour of (f), so the contour must be deformed to 
avoid it, as shown by the solid line. Here p x = 0.001. 
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saddle points. Here the control space is two-dimensional, 
which allows for nontrivial closed loops inside it, and this 
defeats the possibility of attaching any type of label to in¬ 
dividual roots of (37). 

Finally, an interesting consequence of the mixing be¬ 
tween roots is that, at certain specific values of p x and p z , 
the roots must merge, giving double roots of (37). How¬ 
ever, this behaviour depends very sensitively on the mo¬ 
mentum, and it can safely be ignored. In fact, the very dif¬ 
ficulty of tagging the roots, caused by the mixing, makes 
finding the merge momentum an elusive numerical prob¬ 
lem. 


C. Navigating the branch cuts 

Apart from its noticeable effects on the ionization ampli¬ 
tude, the topological change in the temporal branch cuts 
has an obvious effect on the integration contour required to 
produce a correct yield. The closest-approacli times are a 
useful tool in identifying and crossing the branch-cut gates 
when they are present. However, some gates lead to dead¬ 
end regions which do not need to be crossed, which means 
that not all t CA s are useful stepping points. 

Thus, while the correct contour can always be chosen 
by hand given a specific branch cut configuration, it is 
also desirable to have an algorithm that will specify which 
‘slalom gates’ the contour should cross, and in which order 
it should do so. Without such an algorithm, it is impossi¬ 
ble to automate the choice of contour and the calculation 
of photoelectron spectra is impractical. 

Such an algorithm is indeed available, and it is based on 
a geometrical fact shown in Figs. 12(c) and (d): the topo¬ 
logical change in the branch cut connections that happens 
at a soft recollision comes together with a change in the 
sign of the real part of the squared velocity, v(t) 2 , for the 
outer saddle points. 

Thus, in the closed topology of Fig. 12(b), where the 
contour should pass through tc A and . the real kinetic 
energy Re(v(f OA ) 2 ) is positive at those saddle points. In 
contrast, for the open topology of Fig. 12(a), Re(v(tc A ) 2 ) 
and Re(v(fcA ) 2 ) are both negative, and the contour should 
ignore both of those saddle points. 

The physical content of this criterion is quite clear: in the 
quantum orbit formalism, the classically forbidden regions 
are readily identified in the complex time plane as those 
regions where the kinetic energy ^v 2 (t) is negative, or at 
least has negative real part. The undesirable saddle points 
of Fig. 12(a) therefore require the trajectory to tunnel to¬ 
wards the core to be reached. On the other hand, a formal 
proof of the simultaneity of the topological change with the 
emergence of the t CA from the ‘barrier’ is still lacking. 

We display in Fig. 13 some sample integration contours 
produced with this criterion. In general, the navigation is 
straightforward, and there are never any problems when 
p z < 0 or p± is sizeable, in which case contour looks as 
in Fig. 13(a). For near-recolliding electrons, on the other 
hand, the branch cuts do require more careful navigation, 
as we have seen. 

In general, it suffices to take, in order of increasing 
Re(t CA ), those closest-approach times which (i) occur af¬ 
ter ionization, (ii) have reasonably bounded imaginary 




-0.5 -^. . . 1 ■ » 1 . . . . . 

0 71/4 71/2 377/4 77 577/4 377/2 777/4 277 977/4 




FIG. 13. Sample integration contours produced by the f CA 
choosing algorithm described in the text. Most momenta are 
straightforward (a), but near-recollision momenta, like the one 
shown in Fig. 5, do require careful handling, as shown in (b). 
The algorithm correctly handles soft recollisions, (c), as well as 
higher momenta with harder recollisions (d). 

part, and (iii) have positive kinetic energy. To this we 
add two exceptions: the first inward turning point, in 
—7 t/2 < wt < 7t/ 2, which helps avoid crossing regions 
where Re(v(t) 2 ) < 0 where this is not necessary; and the 
first closest-approach time, in 7 t/2 < uit < 37r/2, when it 
can be consistently identified, which can sometimes have 
Re (v(t CA ) 2 ) < 0 at extremely large momenta, but is nev¬ 
ertheless the correct choice in those cases. 

Thus, a concrete set of rules for choosing a contour is 
taking those t CA s for which 

Re(wt CA ) > ut 0 + 7r/5 and 
-|t t < Im(t CA ) < Im(f K ) and 
Re (v(t CA ) 2 ) > -u, 
or — 7t/2 < Re(u/t CA ) < 7t/2 and 
0 < Im(t CA ) < r T , 
or 7r/2 < R e(ojt CA ) < 3 tt/2 and 

Im(t CA ) > 0, 

where u is an adjustable numerical precision, for additional 
flexibility, set by default to 10 -8 a.u. These rules are rela¬ 
tively heuristic and they have a certain amount of leeway 
around them, but they work well over the relevant region 
of photoelectron momentum to produce correct integration 
path choices. 
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Finally, we note that in certain very specific cases, close 
to a soft recollision, the integration path is topologically 
correct but may pass near a singularity of the Coulomb 
kernel; future work aims at avoiding this problem. 


V. RESULTS 

We are now in a position to produce photoelectron spec¬ 
tra. The yield is given by the amplitude in (29), taking 
the integration of (32) to be over the t CA chosen using the 
hopping algorithm described above. The correctness of the 
integration contour is reflected by the lack of discontinuous 
changes in the final integrand, which are easily detected by 
numerical integration routines when they are present. 

In general, the inclusion of the Coulomb correction factor 
(32) has the main effect of increasing the electron yield by 
two or three orders of magnitude compared to the bare SFA 
yield, which is a well known consequence of the Coulomb 
interaction [15,17,19], and is primarily due to the interac¬ 
tions inside the classically-forbidden region. 

This enhancement is particularly strong, and has a par¬ 
ticularly strong variation, near the classical soft recollisions 
p s z . As discussed above, these are accompanied by an ex¬ 
treme enhancement for momenta on one side of p s z and a 
deep and narrow dip on the other. The presence of mul¬ 
tiple saddle points in close proximity, which marks the in¬ 
creased time the electron spends near the core, implies the 
Coulomb correction must be large, but its precise contri¬ 
bution depends sensitively on the details of the topological 
transition at the soft recollision. 

In the neighbourhood of such a soft recollision, the domi¬ 
nating feature is a sharp ridge for very small transverse mo¬ 
menta, as shown in Fig. 14, which terminates at p s z . This 
ridge feature is similar to the cusps detected at low trans¬ 
verse momenta in careful measurements of photoionization 
spectra at long wavelengths [32], which rise suddenly from 
a gaussian background as p± varies, though the measured 
spectra have sharper cusps. 

On the other side of the soft recollisions, i.e. for p z > p s z , 



FIG. 14. Incoherent addition of the sub-cycle ionization yields, 
| (\a(p x ,0,p z )\ 2 + \a(p x ,0, -p 2 )| 2 ), ignoring the shape factor 
R( p), for long wavelength fields (A = 3.1pm, with 7 = 0.31). 
The Coulomb-correction has been integrated over 2.75 laser pe¬ 
riods. Recent experimental results at this wavelength [32] show 
a similar (if more pronounced) peak structure localized at small 
longitudinal momentum. 


this behaviour reverses, and the enhancement turns into a 
strong suppression of ionization. This suppression is ev¬ 
ident in Fig. 14, which contains two transitions closely 
spaced at p s z = 0.063 a.u. and p s z = 0.031 a.u.; the in¬ 
teractions between them create complex structures in their 
immediate vicinity. These fine structures will most likely 
be washed out in any experimental situation, but the spike 
would be expected to remain. 

The rapid change from enhancement to suppression of 
the ionization signal at the topological transitions of Fig. 12 
can be used to show their direct link with the classical 
soft recollisions of Section I, by comparing their respective 
wavelength scaling. To this end we plot in Fig. 15 the 
variation of the on-axis ionization yield - at a small but 
reasonable transverse momentum - with the laser wave¬ 
length. The resulting ridges in the ionization yield closely 
follow the scaling of the classical soft-recollision momenta 
of Fig. 2, shown as red dots, which shows a strong link 
between both structures. Similarly, the opposite scaling of 
the two types of low-energy structure is amenable to ex¬ 
perimental testing. 


VI. CONCLUSIONS 

We have seen how the analytical J?-matrix theory gives 
rise, from first principles, to a quantum-orbit picture of 
strong field ionization which contains the Coulomb inter¬ 
action with the ion as an added action in the phase. In this 
picture the electron trajectory starts at the ion with a real 
position, and as it travels through the classically forbidden 
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FIG. 15. Variation of the on-axis ionization amplitude |a(p)| 2 , 
in an arbitrary logarithmic scale, as a function of the wavelength 
and the corresponding Keldysh parameter 7 . The sudden drops 
in amplitude of Fig. 14 shift along the momentum axis with a 
scaling that closely matches the classical soft-recollision trajec¬ 
tories, shown as red dots. Here the transverse momentum p± 
has been chosen so that the transverse coordinate of the classi¬ 
cal trajectory has a small but positive value, x = 1.07^, at the 
first soft recollision at cut ~ 2ir, to avoid the hard singularity of 
the Coulomb kernel. 
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region it acquires an imaginary component. 

This imaginary component is not generally a problem, 
but it can dominate the trajectory near recollisions, in 
which case the Coulomb interaction can develop a discon¬ 
tinuous change in sign. This marks a branch cut which 
must be carefully handled, by treating the Coulomb in¬ 
teraction Z7(r c i(f)) as an analytic function on the complex 
time plane, with branch cuts inherited from the spatial 
U( r) = —1/yfr^, which the temporal integration contour 
must avoid. 

The key tool for avoiding these branch cuts in the com¬ 
plex time plane are the times of closest approach, which 
satisfy the complex equation r c i(i CA ) • v(t CA ) = 0, and 
which are always present in the middle of each branch-cut 
gate. These times of closest approach have a rich geometry 
of their own, in the complex quantum domain as well as 
in the real-valued simple-man’s model. Moreover, the soft 
recollisions responsible for the Low Energy Structures are 
embedded in this geometry and arise naturally within both 
formalisms. 

Using the times of closest approach, we have developed 
a consistent algorithm which enables the programmatic 
choice of a correct integration contour for all momenta, 
despite topological changes which can be quite brusque as 
the momentum crosses a soft recollision. 

In particular, this formalism incorporates in a natural 
way the known series of trajectories which gives rise to the 
Low-Energy Structures, but it also predicts a second series 
of trajectories at much lower momentum. These trajec¬ 
tories do not appear in theories which neglect the tunnel 
entrance, which explains why they have so far been over¬ 
looked, but they should contribute to the recently discov¬ 
ered Zero-Energy Structures [31-33]. 

More broadly speaking, our work serves as a roadmap 
for the difficulties that must arise in a full first-principles 
semiclassical theory of ionization, and how they might be 
resolved. ARM theory, in its current form, uses only the 
laser-driven trajectory and does not incorporate the ef¬ 


fect of the ion’s Coulomb potential on the trajectory itself. 
However, it is also clear that any first-principles attempt 
to include the Coulomb interaction at that level will be 
subject to the same imaginary-position properties as ARM 
theory, which makes the calculation of the full trajectory 
a nontrivial problem. In particular, such a theory will be 
subject to the same temporal branch cuts as ARM the¬ 
ory, without the benefit of the easily-found ARM closest- 
approach times to help navigate them. 

In addition, by identifying mechanisms which should 
contribute to the observed Zero Energy Structures, we are 
able to propose experiments which can increase its energy 
scale and bring its details within reach of current exper¬ 
imental capabilities. In particular, the energy scaling of 
odd-order soft-recollision trajectories as 7 1 2 3 4 5 6 7 8 I p suggests that 
experiments with harder targets, such as He + , should help 
elucidate the origin of these structures. 

Furthermore, this trajectory-based explanation opens 
the door to testing via schemes which directly alter the 
shape of the trajectory. This includes variation of the pulse 
width, which is known to be a critical variable [26] along 
with the carrier-envelope phase, but it opens the door to 
testing using elliptical pulses or small amounts of second 
or third harmonics to shift the positions of the trajectories 
responsible for the structures. 

This work thus provides a toolbox of approaches with 
which to understand the semiclassical dynamics of low- 
energy photoelectrons in the tunnelling regime, which 
should inform other low-energy ionization phenomena [49— 
51] as well as streaking-based experiments [52] in that 
regime. 
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